import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

data = pd.read_excel('2021改完.xlsx')
data = data[data['Lab Status']!='Negative ID']

def toPos(data):
    fitData = []
    x = data['Latitude'].tolist()
    y = data['Longitude'].tolist()
    for i in range(len(x)):
        fitData.append([x[i],y[i]])
    fitData=np.array(fitData)
    return fitData


from sklearn.cluster import KMeans

fitData=toPos(data)
kmeans=KMeans(n_clusters=4)
kmeans.fit(fitData)

unnest=kmeans.cluster_centers_
nestData=data[data['Lab Status']=='Positive ID']
nest=toPos(nestData)

data=data[data['Lab Status']!='Positive ID']

import math
def snorm(x): # 季节的正态分布
    s2=3.034234234
    s=math.sqrt(s2)
    u=8
    f1=1/(math.sqrt(2*math.pi)*s)
    f2=math.exp(-(x-u)**2/(2*s2))
    return f1*f2

class reportPos:
    def __init__(self,lat,lng,date,GlobalID):
        self.lat=lat
        self.lng=lng
        self.date=date
        self.ID=GlobalID

    def setRRel(self,rel):
        self.rrel=rel

    def caluSRel(self):
        return snorm(self.date.month)

    def caluRel(self):
        if self.rrel==0:
            self.rel=0
        else:
            self.rel = 0.4*self.caluSRel()+0.6*self.rrel

    def output(self):
        print(self.ID,',',self.rel)

dataPos=[]
dataLat=data['Latitude'].tolist()
dataLng=data['Longitude'].tolist()
dataDate=data['Submission Date'].tolist()
dataID=data['GlobalID'].tolist()
for i in range(len(dataLat)):
    dataPos.append(reportPos(dataLat[i],dataLng[i],dataDate[i],dataID[i]))

# 找距离unnest和nest30公里内的点
from math import radians, cos, sin, asin, sqrt

def geodistance(lat1,lng1,lat2,lng2):
    lng1, lat1, lng2, lat2 = map(radians, [float(lng1), float(lat1), float(lng2), float(lat2)]) # 经纬度转换成弧度
    dlon=lng2-lng1
    dlat=lat2-lat1
    a=sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    distance=2*asin(sqrt(a))*6371*1000 # 地球平均半径，6371km
    distance=round(distance/1000,3)
    return distance

unnestSubPos=[]
nestSubPos=[]

for i in unnest:
    result=[]
    lat1,lng1=i
    for j in dataPos:
        dist=geodistance(lat1,lng1,j.lat,j.lng)
        if dist<=30:
            result.append(j)
    unnestSubPos.append(result)

for i in nest:
    result = []
    lat1,lng1=i
    for j in dataPos:
        dist=geodistance(lat1,lng1,j.lat,j.lng)
        if dist<=30:
            result.append(j)
    nestSubPos.append(result)

# 计算正态分布参数
def caluDistr(block,center):
    allLat = [i.lat for i in block]
    allLat.append(center[0])
    allLng = [i.lng for i in block]
    allLng.append(center[1])
    allLat = np.array(allLat)
    allLng = np.array(allLng)
    s1 = allLat.std()**2
    s2 = allLng.std()**2
    u1 = center[0]
    u2 = center[1]
    print(s1,',',s2,',',u1,',',u2)
    return (s1, s2, u1, u2)
    # plt.scatter(allLat, allLng, alpha=0.3)
    # plt.scatter(center[0], center[1], marker='x', color='b')

unnestDistr=[]
nestDistr=[]

for i in range(len(unnestSubPos)):
    block=unnestSubPos[i]
    center=unnest[i]
    unnestDistr.append(caluDistr(block,center))

for i in range(len(nestSubPos)):
    block=nestSubPos[i]
    center=nest[i]
    nestDistr.append(caluDistr(block,center))

# 算每个未标记点的可信度
def getrnorm(s1,s2,u1,u2):
    def r(x,y):
        f1=2*math.pi*u1*u2
        f1=1/f1
        f2=(-1/2) * (((x-u1)**2/s1) + ((y-u2)**2/s2))
        f2=math.exp(f2)
        return f1*f2
    return r

def findNest(j): # 找这个报告属于哪个蜂巢，直接算出概率密度
    minNestRel=lambda x,y:0 # 距离最短蜂巢对应的概率密度函数，如果不更新信度就是0
    minDist=30 # 小于等于30才更新
    for i in range(len(unnest)):
        lat1, lng1 = unnest[i]
        dist=geodistance(lat1,lng1,j.lat,j.lng)
        if dist<=minDist:
            minDist=dist
            distr=unnestDistr[i]
            minNestRel=getrnorm(distr[0],distr[1],distr[2],distr[3])

    for i in range(len(nest)):
        lat1, lng1 = nest[i]
        dist=geodistance(lat1,lng1,j.lat,j.lng)
        if dist<=minDist:
            minDist=dist
            distr=nestDistr[i]
            minNestRel = getrnorm(distr[0], distr[1], distr[2], distr[3])

    j.setRRel(minNestRel(j.lat,j.lng))


for j in dataPos:
    findNest(j)
    j.caluRel()
    j.output()